#This cell sets up basic plotting functions awe
#we will use to visualize the gradient descent routines.

#Make plots interactive
#%matplotlib notebook

#Make plots static
%matplotlib inline

#Make 3D plots
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
#from matplotlib import animation
from IPython.display import HTML
from matplotlib.colors import LogNorm
#from itertools import zip_longest

#Import Numpy
import numpy as np

#Define function for plotting 

def plot_surface(x, y, z, azim=-60, elev=40, dist=10, cmap="RdYlBu_r"):

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    plot_args = {'rstride': 1, 'cstride': 1, 'cmap':cmap,
             'linewidth': 20, 'antialiased': True,
             'vmin': -2, 'vmax': 2}
    ax.plot_surface(x, y, z, **plot_args)
    ax.view_init(azim=azim, elev=elev)
    ax.dist=dist
    ax.set_xlim(-1, 1)
    ax.set_ylim(-1, 1)
    ax.set_zlim(-2, 2)
    
    plt.xticks([-1, -0.5, 0, 0.5, 1], ["-1", "-1/2", "0", "1/2", "1"])
    plt.yticks([-1, -0.5, 0, 0.5, 1], ["-1", "-1/2", "0", "1/2", "1"])
    ax.set_zticks([-2, -1, 0, 1, 2])
    ax.set_zticklabels(["-2", "-1", "0", "1", "2"])
    
    ax.set_xlabel("x", fontsize=18)
    ax.set_ylabel("y", fontsize=18)
    ax.set_zlabel("z", fontsize=18)
    return fig, ax;

def overlay_trajectory_quiver(ax,obj_func,trajectory, color='k'):
    xs=trajectory[:,0]
    ys=trajectory[:,1]
    zs=obj_func(xs,ys)
    ax.quiver(xs[:-1], ys[:-1], zs[:-1], xs[1:]-xs[:-1], ys[1:]-ys[:-1],zs[1:]-zs[:-1],color=color,arrow_length_ratio=0.3)
    
    return ax;

def overlay_trajectory(ax,obj_func,trajectory,label,color='k'):
    xs=trajectory[:,0]
    ys=trajectory[:,1]
    zs=obj_func(xs,ys)
    ax.plot(xs,ys,zs, color, label=label)
    
    return ax;

    
def overlay_trajectory_contour_M(ax,trajectory, label,color='k',lw=2):
    xs=trajectory[:,0]
    ys=trajectory[:,1]
    ax.plot(xs,ys, color, label=label,lw=lw)
    ax.plot(xs[-1],ys[-1],color+'>', markersize=14)
    return ax;

def overlay_trajectory_contour(ax,trajectory, label,color='k',lw=2):
    xs=trajectory[:,0]
    ys=trajectory[:,1]
    ax.plot(xs,ys, color, label=label,lw=lw)
    return ax;
#DEFINE SURFACES WE WILL WORK WITH

#Define monkey saddle and gradient

def saddle(x,y):
    return x**2 - y**2

def grad_saddle(params):
    x=params[0]
    y=params[1]
    grad_x= 2*x
    grad_y= -2*y
    return [grad_x,grad_y]

def monkey_saddle(x,y):
    return x**3 - 3*x*y**2

def grad_monkey_saddle(params):
    x=params[0]
    y=params[1]
    grad_x= 3*x**2-3*y**2
    grad_y= -6*x*y
    return [grad_x,grad_y]

#Define saddle surface

def saddle_surface(x,y,a=1,b=1):
    return a*x**2-b*y**2

def grad_saddle_surface(params,a=1,b=1):
    x=params[0]
    y=params[1]
    grad_x= a*x
    grad_y= -1*b*y
    return [grad_x,grad_y]


# Define minima_surface

def minima_surface(x,y,a=1,b=1):
    return a*x**2+b*y**2-1

def grad_minima_surface(params,a=1,b=1):
    x=params[0]
    y=params[1]
    grad_x= 2*a*x
    grad_y= 2*b*y
    return [grad_x,grad_y]


def beales_function(x,y):
    return np.square(1.5-x+x*y)+np.square(2.25-x+x*y*y)+np.square(2.625-x+x*y**3)
    return f

def grad_beales_function(params):
    x=params[0]
    y=params[1]
    grad_x=2*(1.5-x+x*y)*(-1+y)+2*(2.25-x+x*y**2)*(-1+y**2)+2*(2.625-x+x*y**3)*(-1+y**3)
    grad_y=2*(1.5-x+x*y)*x+4*(2.25-x+x*y**2)*x*y+6*(2.625-x+x*y**3)*x*y**2
    return [grad_x,grad_y]




def contour_beales_function():
    #plot beales function
    x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
    fig, ax = plt.subplots(figsize=(10, 6))
    z=beales_function(x,y)
    cax = ax.contour(x, y, z, levels=np.logspace(0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")
    ax.plot(3,0.5, 'r*', markersize=18)

    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

    ax.set_xlim((-4.5, 4.5))
    ax.set_ylim((-4.5, 4.5))
    
    return fig,ax

def himmel(x,y):
    return( (x**2+y-11)**2 + (x+y**2-7)**2)


def grad_himmel(params):
    x = params[0]
    y = params[1]

    grad_x = 4*x*(x**2+y-11) + 2*(x+y**2-7)
    grad_y = 2*(x**2+y-11) + 4*y*(x+y**2-7)
    return([grad_x,grad_y])





def rosenbrock(x,y):
    a = 1
    b = 100
    return((a-x)**2 + b*(y-x**2)**2)

def grad_rosenbrock(params):
    a = 1
    b = 100
    x=params[0]
    y=params[1]
    grad_x= -2*(a-x) - 4*b*(y-x**2)*x
    grad_y= 2*b*(y-x**2)
    return [grad_x,grad_y]
    
def contour_rosenbrock():
    x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
    fig, ax = plt.subplots(figsize=(10, 6))
    z=rosenbrock(x,y)
    cax = ax.contour(x, y, z, levels=np.logspace(0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")
    ax.plot(1,1, 'r*', markersize=18)

    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

    ax.set_xlim((-4.5, 4.5))
    ax.set_ylim((-4.5, 4.5))
    
    return fig,ax

def contour_himmel():
    x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
    fig, ax = plt.subplots(figsize=(10, 6))
    z=himmel(x,y)
    cax = ax.contour(x, y, z, levels=np.logspace(0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")

    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

    ax.set_xlim((-4.5, 4.5))
    ax.set_ylim((-4.5, 4.5))
    
    return fig,ax

#Make plots of surfaces
plt.close() # closes previous plots
x, y = np.mgrid[-1:1:31j, -1:1:31j]
fig1,ax1=plot_surface(x,y,monkey_saddle(x,y))
fig2,ax2=plot_surface(x,y,saddle_surface(x,y))
fig3,ax3=plot_surface(x,y,minima_surface(x,y,5),0)

#Contour plot of Beale's Function

fig4,ax4 =contour_beales_function()
plt.show()

def contour_saddle():
    x, y = np.meshgrid(np.arange(-2.5, 2.5, 0.2), np.arange(-2.5, 2.5, 0.2))
    fig, ax = plt.subplots(figsize=(10, 6))
    z= saddle(x,y)
    cax = ax.contour(x, y, z,levels=np.linspace(-10,10,35), norm=LogNorm(), cmap="RdYlBu_r")

    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

    ax.set_xlim((-2.5, 2.5))
    ax.set_ylim((-2.5, 2.5))
    
    return fig,ax

#This writes a simple gradient descent, gradient descent+ momentum,
#nesterov. 

#Mean-gradient based methods
def gd(grad, init, n_epochs=1000, eta=10**-4, noise_strength=0):
    #This is a simple optimizer
    params=np.array(init)
    param_traj=np.zeros([n_epochs+1,2])
    param_traj[0,]=init
    v=0;
    for j in range(n_epochs):
        noise=noise_strength*np.random.randn(params.size)
        v=eta*(np.array(grad(params))+noise)
        params=params-v
        param_traj[j+1,]=params
    return param_traj


def gd_with_mom(grad, init, n_epochs=5000, eta=10**-4, gamma=0.9,noise_strength=0):
    params=np.array(init)
    param_traj=np.zeros([n_epochs+1,2])
    param_traj[0,]=init
    v=0
    for j in range(n_epochs):
        noise=noise_strength*np.random.randn(params.size)
        v=gamma*v+eta*(np.array(grad(params))+noise)
        params=params-v
        param_traj[j+1,]=params
    return param_traj

def NAG(grad, init, n_epochs=5000, eta=10**-4, gamma=0.9,noise_strength=0):
    params=np.array(init)
    param_traj=np.zeros([n_epochs+1,2])
    param_traj[0,]=init
    v=0
    for j in range(n_epochs):
        noise=noise_strength*np.random.randn(params.size)
        params_nesterov=params-gamma*v
        v=gamma*v+eta*(np.array(grad(params_nesterov))+noise)
        params=params-v
        param_traj[j+1,]=params
    return param_traj
# Investigate effect of learning rate in GD
plt.close()
a,b = 1.0,1.0
x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
fig, ax = plt.subplots(figsize=(10, 6))
z=np.abs(minima_surface(x,y,a,b))
ax.contour(x, y, z, levels=np.logspace(0.0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")
ax.plot(0,0, 'r*', markersize=18)

#initial point
init1=[-2,4]
init2=[-1.7,4]
init3=[-1.5,4]
init4=[-3,4.5]
eta1=0.1
eta2=0.5
eta3=1
eta4=1.01
gd_1=gd(grad_minima_surface,init1, n_epochs=100, eta=eta1)
gd_2=gd(grad_minima_surface,init2, n_epochs=100, eta=eta2)
gd_3=gd(grad_minima_surface,init3, n_epochs=100, eta=eta3)
gd_4=gd(grad_minima_surface,init4, n_epochs=10, eta=eta4)
#print(gd_1)
overlay_trajectory_contour(ax,gd_1,'$\eta=$%s'% eta1,'g--*', lw=0.5)
overlay_trajectory_contour(ax,gd_2,'$\eta=$%s'% eta2,'b-<', lw=0.5)
overlay_trajectory_contour(ax,gd_3,'$\eta=$%s'% eta3,'->', lw=0.5)
overlay_trajectory_contour(ax,gd_4,'$\eta=$%s'% eta4,'c-o', lw=0.5)
plt.legend(loc=2)
plt.show()
fig.savefig("GD3regimes.pdf", bbox_inches='tight')

# Investigate effect of learning rate in GD
plt.close()
a,b = 1.0,1.0
x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
fig, ax = plt.subplots(figsize=(10, 6))
z=np.abs(minima_surface(x,y,a,b))
ax.contour(x, y, z, levels=np.logspace(0.0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")
ax.plot(0,0, 'r*', markersize=18)

#initial point
init1=[-2,4]
init2=[-1.7,4]
init3=[-1.5,4]
init4=[-3,4.5]
eta1=0.1
eta2=0.5
eta3=1
eta4=1.01
gd_1=gd_with_mom(grad_minima_surface,init1, n_epochs=100, eta=eta1)
gd_2=gd_with_mom(grad_minima_surface,init2, n_epochs=100, eta=eta2)
gd_3=gd_with_mom(grad_minima_surface,init3, n_epochs=100, eta=eta3)
gd_4=gd_with_mom(grad_minima_surface,init4, n_epochs=10, eta=eta4)
#print(gd_1)
overlay_trajectory_contour(ax,gd_1,'$\eta=$%s'% eta1,'g--*', lw=0.5)
overlay_trajectory_contour(ax,gd_2,'$\eta=$%s'% eta2,'b-<', lw=0.5)
overlay_trajectory_contour(ax,gd_3,'$\eta=$%s'% eta3,'->', lw=0.5)
overlay_trajectory_contour(ax,gd_4,'$\eta=$%s'% eta4,'c-o', lw=0.5)
plt.legend(loc=2)
plt.show()
#fig.savefig("GD3regimes.pdf", bbox_inches='tight')

# Investigate effect of learning rate in GD
plt.close()
a,b = 1.0,1.0
x, y = np.meshgrid(np.arange(-4.5, 4.5, 0.2), np.arange(-4.5, 4.5, 0.2))
fig, ax = plt.subplots(figsize=(10, 6))
z=np.abs(minima_surface(x,y,a,b))
ax.contour(x, y, z, levels=np.logspace(0.0, 5, 35), norm=LogNorm(), cmap="RdYlBu_r")
ax.plot(0,0, 'r*', markersize=18)

#initial point
init1=[-2,4]
init2=[-1.7,4]
init3=[-1.5,4]
init4=[-3,4.5]
eta1=0.1
eta2=0.5
eta3=1
eta4=1.01
gd_1=NAG(grad_minima_surface,init1, n_epochs=100, eta=eta1)
gd_2=NAG(grad_minima_surface,init2, n_epochs=100, eta=eta2)
#gd_3=NAG(grad_minima_surface,init3, n_epochs=100, eta=eta3)
#gd_4=NAG(grad_minima_surface,init4, n_epochs=10, eta=eta4)
#print(gd_1)
overlay_trajectory_contour(ax,gd_1,'$\eta=$%s'% eta1,'g--*', lw=0.5)
overlay_trajectory_contour(ax,gd_2,'$\eta=$%s'% eta2,'b-<', lw=0.5)
#overlay_trajectory_contour(ax,gd_3,'$\eta=$%s'% eta3,'->', lw=0.5)
#overlay_trajectory_contour(ax,gd_4,'$\eta=$%s'% eta4,'c-o', lw=0.5)
plt.legend(loc=2)
plt.show()
#fig.savefig("GD3regimes.pdf", bbox_inches='tight')

################################################################################
# Methods that exploit first and second moments of gradient: RMS-PROP and ADAMS
################################################################################

def rms_prop(grad, init, n_epochs=5000, eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0):
    params=np.array(init)
    param_traj=np.zeros([n_epochs+1,2])
    param_traj[0,]=init#Import relevant packages
    grad_sq=0;
    for j in range(n_epochs):
        noise=noise_strength*np.random.randn(params.size)
        g=np.array(grad(params))+noise
        grad_sq=beta*grad_sq+(1-beta)*g*g
        v=eta*np.divide(g,np.sqrt(grad_sq+epsilon))
        params= params-v
        param_traj[j+1,]=params
    return param_traj
                        
                        
def adams(grad, init, n_epochs=5000, eta=10**-4, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0):
    params=np.array(init)
    param_traj=np.zeros([n_epochs+1,2])
    param_traj[0,]=init
    v=0;
    grad_sq=0;
    for j in range(n_epochs):
        noise=noise_strength*np.random.randn(params.size)
        g=np.array(grad(params))+noise
        v=gamma*v+(1-gamma)*g
        grad_sq=beta*grad_sq+(1-beta)*g*g
        v_hat=v/(1-gamma**(j+1))
        grad_sq_hat=grad_sq/(1-beta**(j+1))
        params=params-eta*np.divide(v_hat,np.sqrt(grad_sq_hat+epsilon))
        param_traj[j+1,]=params
    return param_traj
plt.close()
#Make static plot of the results
Nsteps=10**4
lr_l=10**-3
lr_s=10**-6

init1=np.array([4,3])
fig1, ax1=contour_beales_function()

gd_trajectory1=gd(grad_beales_function,init1,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory1=gd_with_mom(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory1=NAG(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory1=rms_prop(grad_beales_function,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory1=adams(grad_beales_function,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')

init3=np.array([-1,4])

gd_trajectory3=gd(grad_beales_function,init3,10**5, eta=lr_s, noise_strength=0)
gdm_trajectory3=gd_with_mom(grad_beales_function,init3,10**5,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory3=NAG(grad_beales_function,init3,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory3=rms_prop(grad_beales_function,init3,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory3=adams(grad_beales_function,init3,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory3, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory3, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory3, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory3,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory3,'ADAMS', 'r')

init4=np.array([-2,-4])

gd_trajectory4=gd(grad_beales_function,init4,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory4=gd_with_mom(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory4=NAG(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory4=rms_prop(grad_beales_function,init4,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory4=adams(grad_beales_function,init4,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory4, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory4, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory4, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory4,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory4,'ADAMS', 'r')

plt.show()

plt.close()
#Make static plot of the results
Nsteps=10**4
lr_l=10**-6
lr_s=10**-6

init1=np.array([4,3])
fig1, ax1=contour_beales_function()

gd_trajectory1=gd(grad_beales_function,init1,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory1=gd_with_mom(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory1=NAG(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory1=rms_prop(grad_beales_function,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory1=adams(grad_beales_function,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')

init3=np.array([-1,4])

gd_trajectory3=gd(grad_beales_function,init3,10**5, eta=lr_s, noise_strength=0)
gdm_trajectory3=gd_with_mom(grad_beales_function,init3,10**5,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory3=NAG(grad_beales_function,init3,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory3=rms_prop(grad_beales_function,init3,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory3=adams(grad_beales_function,init3,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory3, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory3, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory3, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory3,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory3,'ADAMS', 'r')

init4=np.array([-2,-4])

gd_trajectory4=gd(grad_beales_function,init4,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory4=gd_with_mom(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory4=NAG(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory4=rms_prop(grad_beales_function,init4,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory4=adams(grad_beales_function,init4,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory4, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory4, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory4, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory4,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory4,'ADAMS', 'r')

plt.show()

plt.close()
#Make static plot of the results
Nsteps=10**4
lr_l=10**-3
lr_s=10**-3

init1=np.array([4,3])
fig1, ax1=contour_beales_function()

gd_trajectory1=gd(grad_beales_function,init1,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory1=gd_with_mom(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory1=NAG(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory1=rms_prop(grad_beales_function,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory1=adams(grad_beales_function,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')

init3=np.array([-1,4])

gd_trajectory3=gd(grad_beales_function,init3,10**5, eta=lr_s, noise_strength=0)
gdm_trajectory3=gd_with_mom(grad_beales_function,init3,10**5,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory3=NAG(grad_beales_function,init3,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory3=rms_prop(grad_beales_function,init3,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory3=adams(grad_beales_function,init3,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory3, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory3, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory3, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory3,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory3,'ADAMS', 'r')

init4=np.array([-2,-4])

gd_trajectory4=gd(grad_beales_function,init4,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory4=gd_with_mom(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory4=NAG(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory4=rms_prop(grad_beales_function,init4,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=0)
adam_trajectory4=adams(grad_beales_function,init4,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)

overlay_trajectory_contour_M(ax1,gd_trajectory4, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory4, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory4, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory4,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory4,'ADAMS', 'r')

plt.show()
/tmp/ipykernel_739963/3449801150.py:47: RuntimeWarning: overflow encountered in scalar multiply
  grad_x=2*(1.5-x+x*y)*(-1+y)+2*(2.25-x+x*y**2)*(-1+y**2)+2*(2.625-x+x*y**3)*(-1+y**3)
/tmp/ipykernel_739963/3449801150.py:47: RuntimeWarning: overflow encountered in scalar power
  grad_x=2*(1.5-x+x*y)*(-1+y)+2*(2.25-x+x*y**2)*(-1+y**2)+2*(2.625-x+x*y**3)*(-1+y**3)
/tmp/ipykernel_739963/3449801150.py:48: RuntimeWarning: overflow encountered in scalar multiply
  grad_y=2*(1.5-x+x*y)*x+4*(2.25-x+x*y**2)*x*y+6*(2.625-x+x*y**3)*x*y**2
/tmp/ipykernel_739963/3449801150.py:48: RuntimeWarning: overflow encountered in scalar power
  grad_y=2*(1.5-x+x*y)*x+4*(2.25-x+x*y**2)*x*y+6*(2.625-x+x*y**3)*x*y**2
/tmp/ipykernel_739963/3449801150.py:47: RuntimeWarning: invalid value encountered in scalar add
  grad_x=2*(1.5-x+x*y)*(-1+y)+2*(2.25-x+x*y**2)*(-1+y**2)+2*(2.625-x+x*y**3)*(-1+y**3)
/tmp/ipykernel_739963/3449801150.py:48: RuntimeWarning: invalid value encountered in scalar add
  grad_y=2*(1.5-x+x*y)*x+4*(2.25-x+x*y**2)*x*y+6*(2.625-x+x*y**3)*x*y**2

plt.close()
#Make static plot of the results
Nsteps=10**5
lr_l=10**-3
lr_s=10**-4
noise = 400.0
init1=np.array([4,3])
fig1, ax1=contour_beales_function()

gd_trajectory1=gd(grad_beales_function,init1,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectory1=gd_with_mom(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)
#NAG_trajectory1=NAG(grad_beales_function,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
#rms_prop_trajectory1=rms_prop(grad_beales_function,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
#adam_trajectory1=adams(grad_beales_function,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')

init3=np.array([-1,4])

gd_trajectory3=gd(grad_beales_function,init3,10**5, eta=lr_s, noise_strength=noise)
gdm_trajectory3=gd_with_mom(grad_beales_function,init3,10**5,eta=lr_s, gamma=0.9,noise_strength=noise)
#NAG_trajectory3=NAG(grad_beales_function,init3,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
#rms_prop_trajectory3=rms_prop(grad_beales_function,init3,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
#adam_trajectory3=adams(grad_beales_function,init3,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory3, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory3, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory3, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory3,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory3,'ADAMS', 'r')

init4=np.array([-2,-4])

gd_trajectory4=gd(grad_beales_function,init4,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectory4=gd_with_mom(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)
#NAG_trajectory4=NAG(grad_beales_function,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
#rms_prop_trajectory4=rms_prop(grad_beales_function,init4,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
#adam_trajectory4=adams(grad_beales_function,init4,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory4, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory4, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory4, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory4,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory4,'ADAMS', 'r')

plt.show()


#Make plots of surfaces
plt.close() # closes previous plots
x, y = np.mgrid[-1:1:31j, -1:1:31j]
fig1,ax1=contour_rosenbrock()
plt.show()

plt.close()
#Make static plot of the results
Nsteps=10**5
lr_l=10**-3
lr_s=10**-3
noise = 0.0
init1=np.array([4,3])
fig1, ax1=contour_rosenbrock()

gd_trajectory1=gd(grad_rosenbrock,init1,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectory1=gd_with_mom(grad_rosenbrock,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)
NAG_trajectory1=NAG(grad_rosenbrock,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory1=rms_prop(grad_rosenbrock,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
adam_trajectory1=adams(grad_rosenbrock,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')

init3=np.array([-1,4])

gd_trajectory3=gd(grad_rosenbrock,init3,10**5, eta=lr_s, noise_strength=noise)
gdm_trajectory3=gd_with_mom(grad_rosenbrock,init3,10**5,eta=lr_s, gamma=0.9,noise_strength=noise)
NAG_trajectory3=NAG(grad_rosenbrock,init3,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory3=rms_prop(grad_rosenbrock,init3,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
adam_trajectory3=adams(grad_rosenbrock,init3,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory3, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory3, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory3, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory3,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory3,'ADAMS', 'r')

init4=np.array([-2,-4])

gd_trajectory4=gd(grad_rosenbrock,init4,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectory4=gd_with_mom(grad_rosenbrock,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)
NAG_trajectory4=NAG(grad_rosenbrock,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory4=rms_prop(grad_rosenbrock,init4,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
adam_trajectory4=adams(grad_rosenbrock,init4,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory4, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory4, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory4, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory4,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory4,'ADAMS', 'r')

plt.show()
/tmp/ipykernel_739963/3753588982.py:96: RuntimeWarning: overflow encountered in scalar power
  grad_x= -2*(a-x) - 4*b*(y-x**2)*x
/tmp/ipykernel_739963/3753588982.py:97: RuntimeWarning: overflow encountered in scalar power
  grad_y= 2*b*(y-x**2)
/tmp/ipykernel_739963/3753588982.py:96: RuntimeWarning: invalid value encountered in scalar subtract
  grad_x= -2*(a-x) - 4*b*(y-x**2)*x
/tmp/ipykernel_739963/3753588982.py:97: RuntimeWarning: invalid value encountered in scalar subtract
  grad_y= 2*b*(y-x**2)
/tmp/ipykernel_739963/3753588982.py:96: RuntimeWarning: overflow encountered in scalar multiply
  grad_x= -2*(a-x) - 4*b*(y-x**2)*x
/tmp/ipykernel_739963/4156107574.py:39: RuntimeWarning: invalid value encountered in add
  v=gamma*v+eta*(np.array(grad(params_nesterov))+noise)
/tmp/ipykernel_739963/4156107574.py:14: RuntimeWarning: invalid value encountered in subtract
  params=params-v

plt.plot(np.log10(np.linalg.norm(gdm_trajectory4 - [1,1],axis=1)))
plt.plot(np.log10(np.linalg.norm(gd_trajectory4 - [1,1],axis=1)))
plt.plot(np.log10(np.linalg.norm(adam_trajectory4 - [1,1],axis=1)))

contour_himmel()

plt.close()
#Make static plot of the results
Nsteps=10**2
lr_l=10**-2
lr_s=10**-4
noise = 0.0
init1=np.array([-3,0])
fig1, ax1=contour_himmel()

gd_trajectory1=gd(grad_himmel,init1,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectory1=gd_with_mom(grad_himmel,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)
NAG_trajectory1=NAG(grad_himmel,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory1=rms_prop(grad_himmel,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
adam_trajectory1=adams(grad_himmel,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')

init3=np.array([-1,4])

gd_trajectory3=gd(grad_himmel,init3,10**5, eta=lr_s, noise_strength=noise)
gdm_trajectory3=gd_with_mom(grad_himmel,init3,10**5,eta=lr_s, gamma=0.9,noise_strength=noise)
NAG_trajectory3=NAG(grad_himmel,init3,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory3=rms_prop(grad_himmel,init3,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
adam_trajectory3=adams(grad_himmel,init3,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory3, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory3, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory3, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory3,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory3,'ADAMS', 'r')

init4=np.array([-2,-4])

gd_trajectory4=gd(grad_himmel,init4,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectory4=gd_with_mom(grad_himmel,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)
NAG_trajectory4=NAG(grad_himmel,init4,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory4=rms_prop(grad_himmel,init4,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
adam_trajectory4=adams(grad_himmel,init4,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory4, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory4, 'GDM','m')
overlay_trajectory_contour_M(ax1,NAG_trajectory4, 'NAG','c--')
overlay_trajectory_contour_M(ax1,rms_prop_trajectory4,'RMS', 'b-.')
overlay_trajectory_contour_M(ax1,adam_trajectory4,'ADAMS', 'r')

plt.show()

plt.plot(np.log10(np.linalg.norm(gdm_trajectory1 - [3,2],axis=1)))
plt.plot(np.log10(np.linalg.norm(gd_trajectory1 - [3,2],axis=1)))
plt.plot(np.log10(np.linalg.norm(adam_trajectory1 - [3,2],axis=1)))
/tmp/ipykernel_739963/3398975203.py:3: RuntimeWarning: divide by zero encountered in log10
  plt.plot(np.log10(np.linalg.norm(adam_trajectory1 - [3,2],axis=1)))

plt.plot(gdm_trajectory1-[gdm_trajectory1[-1,-1],gdm_trajectory1[-1,-1]])

plt.close()
#Make static plot of the results
Nsteps=10**3
lr_l=10**-2
lr_s=10**-4
noise = 10.0
init1=np.array([-3,0])
fig1, ax1=contour_himmel()

gd_trajectory1=gd(grad_himmel,init1,Nsteps, eta=lr_s, noise_strength=0)
gdm_trajectory1=gd_with_mom(grad_himmel,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
NAG_trajectory1=NAG(grad_himmel,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
rms_prop_trajectory1=rms_prop(grad_himmel,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
adam_trajectory1=adams(grad_himmel,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)


gd_trajectoryNoise=gd(grad_himmel,init1,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectoryNoise=gd_with_mom(grad_himmel,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
#overlay_trajectory_contour_M(ax1,gdm_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,gd_trajectoryNoise, 'GD_N','c-')
#overlay_trajectory_contour_M(ax1,gdm_trajectoryNoise, 'GDM_N','r-')
#overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')


plt.show()

plt.close()
#Make static plot of the results
Nsteps=10**5
lr_l=10**-2
lr_s=10**-4
noise = 3.0
init1=np.array([-2,0])
fig1, ax1=contour_saddle()

gd_trajectory1=gd(grad_saddle,init1,Nsteps, eta=lr_s, noise_strength=0.1)
gdm_trajectory1=gd_with_mom(grad_saddle,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
#NAG_trajectory1=NAG(grad_saddle,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=0)
#rms_prop_trajectory1=rms_prop(grad_saddle,init1,Nsteps,eta=lr_l, beta=0.9,epsilon=10**-8,noise_strength=noise)
#adam_trajectory1=adams(grad_saddle,init1,Nsteps,eta=lr_l, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=noise)


gd_trajectoryNoise=gd(grad_saddle,init1,Nsteps, eta=lr_s, noise_strength=noise)
gdm_trajectoryNoise=gd_with_mom(grad_saddle,init1,Nsteps,eta=lr_s, gamma=0.9,noise_strength=noise)

overlay_trajectory_contour_M(ax1,gd_trajectory1, 'GD','k')
overlay_trajectory_contour_M(ax1,gdm_trajectory1, 'GD-M','r')

#overlay_trajectory_contour_M(ax1,gdm_trajectory1, 'GDM','m')
overlay_trajectory_contour_M(ax1,gd_trajectoryNoise, 'GD_N','c-')
#overlay_trajectory_contour_M(ax1,gdm_trajectoryNoise, 'GDM_N','r-')
#overlay_trajectory_contour_M(ax1,NAG_trajectory1, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory1,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory1,'ADAMS', 'r')

plt.legend(loc=2)

#init2=np.array([1.5,1.5])
#gd_trajectory2=gd(grad_beales_function,init2,Nsteps, eta=10**-6, noise_strength=0)
#gdm_trajectory2=gd_with_mom(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#NAG_trajectory2=NAG(grad_beales_function,init2,Nsteps,eta=10**-6, gamma=0.9,noise_strength=0)
#rms_prop_trajectory2=rms_prop(grad_beales_function,init2,Nsteps,eta=10**-3, beta=0.9,epsilon=10**-8,noise_strength=0)
#adam_trajectory2=adams(grad_beales_function,init2,Nsteps,eta=10**-3, gamma=0.9, beta=0.99,epsilon=10**-8,noise_strength=0)
#overlay_trajectory_contour_M(ax1,gdm_trajectory2, 'GDM','m')
#overlay_trajectory_contour_M(ax1,NAG_trajectory2, 'NAG','c--')
#overlay_trajectory_contour_M(ax1,rms_prop_trajectory2,'RMS', 'b-.')
#overlay_trajectory_contour_M(ax1,adam_trajectory2,'ADAMS', 'r')


plt.show()

def contour_saddle():
    x, y = np.meshgrid(np.arange(-2.5, 2.5, 0.2), np.arange(-2.5, 2.5, 0.2))
    fig, ax = plt.subplots(figsize=(10, 6))
    z= saddle(x,y)
    cax = ax.contourf(x, y, z,levels=np.linspace(-10,10,35), cmap="bwr")

    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')

    ax.set_xlim((-2.5, 2.5))
    ax.set_ylim((-2.5, 2.5))
    
    return fig,ax
contour_saddle()

saddle(2,0)
4

x, y = np.mgrid[-1:1:31j, -1:1:31j]
fig1,ax1=plot_surface(x,y,saddle(x,y))

def plot_surface(x, y, z, azim=0, elev=60, dist=30, cmap="RdYlBu_r"):

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    plot_args = {'rstride': 1, 'cstride': 1, 'cmap':cmap,
             'linewidth': 20, 'antialiased': True,
             'vmin': 0, 'vmax': 200}
    ax.plot_surface(x, y, z, **plot_args)
    ax.view_init(azim=azim, elev=elev)
    ax.dist=dist
    ax.set_xlim(-1, 1)
    ax.set_ylim(-1, 1)
    ax.set_zlim(-2, 2)
    
    plt.xticks([-5, -2, 0, 2, 5], ["-1", "-1/2", "0", "1/2", "1"])
    plt.yticks([-5, -2, 0, 2, 5], ["-1", "-1/2", "0", "1/2", "1"])
    ax.set_zticks([0, 50, 100, 150,200])
    ax.set_zticklabels(["0", "50", "100", "150","200"])
    
    ax.set_xlabel("x", fontsize=18)
    ax.set_ylabel("y", fontsize=18)
    ax.set_zlabel("z", fontsize=18)
    return fig, ax;

himmel(x,y)
array([[2.50000000e+02, 1.82716049e+02, 1.39382716e+02, 1.16000000e+02,
        1.08864198e+02, 1.14567901e+02, 1.30000000e+02, 1.52345679e+02,
        1.79086420e+02, 2.08000000e+02, 2.37160494e+02, 2.64938272e+02,
        2.90000000e+02, 3.11308642e+02, 3.28123457e+02, 3.40000000e+02,
        3.46790123e+02, 3.48641975e+02, 3.46000000e+02, 3.39604938e+02,
        3.30493827e+02, 3.20000000e+02, 3.09753086e+02, 3.01679012e+02,
        2.98000000e+02, 3.01234568e+02, 3.14197531e+02, 3.40000000e+02,
        3.82049383e+02, 4.44049383e+02, 5.30000000e+02],
       [2.11160494e+02, 1.39580247e+02, 9.20987654e+01, 6.47160494e+01,
        5.37283951e+01, 5.57283951e+01, 6.76049383e+01, 8.65432099e+01,
        1.10024691e+02, 1.35827160e+02, 1.62024691e+02, 1.86987654e+02,
        2.09382716e+02, 2.28172840e+02, 2.42617284e+02, 2.52271605e+02,
        2.56987654e+02, 2.56913580e+02, 2.52493827e+02, 2.44469136e+02,
        2.33876543e+02, 2.22049383e+02, 2.10617284e+02, 2.01506173e+02,
        1.96938272e+02, 1.99432099e+02, 2.11802469e+02, 2.37160494e+02,
        2.78913580e+02, 3.40765432e+02, 4.26716049e+02],
       [1.94493827e+02, 1.18765432e+02, 6.72839506e+01, 3.60493827e+01,
        2.13580247e+01, 1.98024691e+01, 2.82716049e+01, 4.39506173e+01,
        6.43209877e+01, 8.71604938e+01, 1.10543210e+02, 1.32839506e+02,
        1.52716049e+02, 1.69135802e+02, 1.81358025e+02, 1.88938272e+02,
        1.91728395e+02, 1.89876543e+02, 1.83827160e+02, 1.74320988e+02,
        1.62395062e+02, 1.49382716e+02, 1.36913580e+02, 1.26913580e+02,
        1.21604938e+02, 1.23506173e+02, 1.35432099e+02, 1.60493827e+02,
        2.02098765e+02, 2.63950617e+02, 3.50049383e+02],
       [1.96000000e+02, 1.16271605e+02, 6.09382716e+01, 2.60000000e+01,
        7.75308642e+00, 2.79012346e+00, 8.00000000e+00, 2.05679012e+01,
        3.79753086e+01, 5.80000000e+01, 7.87160494e+01, 9.84938272e+01,
        1.16000000e+02, 1.30197531e+02, 1.40345679e+02, 1.46000000e+02,
        1.47012346e+02, 1.43530864e+02, 1.36000000e+02, 1.25160494e+02,
        1.12049383e+02, 9.80000000e+01, 8.46419753e+01, 7.39012346e+01,
        6.80000000e+01, 6.94567901e+01, 8.10864198e+01, 1.06000000e+02,
        1.47604938e+02, 2.09604938e+02, 2.96000000e+02],
       [2.11975309e+02, 1.28395062e+02, 6.93580247e+01, 3.08641975e+01,
        9.20987654e+00, 9.87654321e-01, 3.08641975e+00, 1.26913580e+01,
        2.72839506e+01, 4.46419753e+01, 6.28395062e+01, 8.02469136e+01,
        9.55308642e+01, 1.07654321e+02, 1.15876543e+02, 1.19753086e+02,
        1.19135802e+02, 1.14172840e+02, 1.05308642e+02, 9.32839506e+01,
        7.91358025e+01, 6.41975309e+01, 5.00987654e+01, 3.87654321e+01,
        3.24197531e+01, 3.35802469e+01, 4.50617284e+01, 6.99753086e+01,
        1.11728395e+02, 1.74024691e+02, 2.60864198e+02],
       [2.39012346e+02, 1.51728395e+02, 8.91358025e+01, 4.72345679e+01,
        2.23209877e+01, 1.09876543e+01, 1.01234568e+01, 1.69135802e+01,
        2.88395062e+01, 4.36790123e+01, 5.95061728e+01, 7.46913580e+01,
        8.79012346e+01, 9.80987654e+01, 1.04543210e+02, 1.06790123e+02,
        1.04691358e+02, 9.83950617e+01, 8.83456790e+01, 7.52839506e+01,
        6.02469136e+01, 4.45679012e+01, 2.98765432e+01, 1.80987654e+01,
        1.14567901e+01, 1.24691358e+01, 2.39506173e+01, 4.90123457e+01,
        9.10617284e+01, 1.53802469e+02, 2.41234568e+02],
       [2.74000000e+02, 1.83160494e+02, 1.17160494e+02, 7.20000000e+01,
        4.39753086e+01, 2.96790123e+01, 2.60000000e+01, 3.01234568e+01,
        3.95308642e+01, 5.20000000e+01, 6.56049383e+01, 7.87160494e+01,
        9.00000000e+01, 9.84197531e+01, 1.03234568e+02, 1.04000000e+02,
        1.00567901e+02, 9.30864198e+01, 8.20000000e+01, 6.80493827e+01,
        5.22716049e+01, 3.60000000e+01, 2.08641975e+01, 8.79012346e+00,
        2.00000000e+00, 3.01234568e+00, 1.46419753e+01, 4.00000000e+01,
        8.24938272e+01, 1.45827160e+02, 2.34000000e+02],
       [3.14123457e+02, 2.19876543e+02, 1.50617284e+02, 1.02345679e+02,
        7.13580247e+01, 5.42469136e+01, 4.79012346e+01, 4.95061728e+01,
        5.65432099e+01, 6.67901235e+01, 7.83209877e+01, 8.95061728e+01,
        9.90123457e+01, 1.05802469e+02, 1.09135802e+02, 1.08567901e+02,
        1.03950617e+02, 9.54320988e+01, 8.34567901e+01, 6.87654321e+01,
        5.23950617e+01, 3.56790123e+01, 2.02469136e+01, 8.02469136e+00,
        1.23456790e+00, 2.39506173e+00, 1.43209877e+01, 4.01234568e+01,
        8.32098765e+01, 1.47283951e+02, 2.36345679e+02],
       [3.56864198e+02, 2.59358025e+02, 1.86987654e+02, 1.35753086e+02,
        1.01950617e+02, 8.21728395e+01, 7.33086420e+01, 7.25432099e+01,
        7.73580247e+01, 8.55308642e+01, 9.51358025e+01, 1.04543210e+02,
        1.12419753e+02, 1.17728395e+02, 1.19728395e+02, 1.17975309e+02,
        1.12320988e+02, 1.02913580e+02, 9.01975309e+01, 7.49135802e+01,
        5.80987654e+01, 4.10864198e+01, 2.55061728e+01, 1.32839506e+01,
        6.64197531e+00, 8.09876543e+00, 2.04691358e+01, 4.68641975e+01,
        9.06913580e+01, 1.55654321e+02, 2.45753086e+02],
       [4.00000000e+02, 2.99382716e+02, 2.24049383e+02, 1.70000000e+02,
        1.33530864e+02, 1.11234568e+02, 1.00000000e+02, 9.70123457e+01,
        9.97530864e+01, 1.06000000e+02, 1.13827160e+02, 1.21604938e+02,
        1.28000000e+02, 1.31975309e+02, 1.32790123e+02, 1.30000000e+02,
        1.23456790e+02, 1.13308642e+02, 1.00000000e+02, 8.42716049e+01,
        6.71604938e+01, 5.00000000e+01, 3.44197531e+01, 2.23456790e+01,
        1.60000000e+01, 1.79012346e+01, 3.08641975e+01, 5.80000000e+01,
        1.02716049e+02, 1.68716049e+02, 2.60000000e+02],
       [4.41604938e+02, 3.38024691e+02, 2.59876543e+02, 2.03160494e+02,
        1.64172840e+02, 1.39506173e+02, 1.26049383e+02, 1.20987654e+02,
        1.21802469e+02, 1.26271605e+02, 1.32469136e+02, 1.38765432e+02,
        1.43827160e+02, 1.46617284e+02, 1.46395062e+02, 1.42716049e+02,
        1.35432099e+02, 1.24691358e+02, 1.10938272e+02, 9.49135802e+01,
        7.76543210e+01, 6.04938272e+01, 4.50617284e+01, 3.32839506e+01,
        2.73827160e+01, 2.98765432e+01, 4.35802469e+01, 7.16049383e+01,
        1.17358025e+02, 1.84543210e+02, 2.77160494e+02],
       [4.80049383e+02, 3.73654321e+02, 2.92839506e+02, 2.33604938e+02,
        1.92246914e+02, 1.65358025e+02, 1.49827160e+02, 1.42839506e+02,
        1.41876543e+02, 1.44716049e+02, 1.49432099e+02, 1.54395062e+02,
        1.58271605e+02, 1.60024691e+02, 1.58913580e+02, 1.54493827e+02,
        1.46617284e+02, 1.35432099e+02, 1.21382716e+02, 1.05209877e+02,
        8.79506173e+01, 7.09382716e+01, 5.58024691e+01, 4.44691358e+01,
        3.91604938e+01, 4.23950617e+01, 5.69876543e+01, 8.60493827e+01,
        1.32987654e+02, 2.01506173e+02, 2.95604938e+02],
       [5.14000000e+02, 4.04938272e+02, 3.21604938e+02, 2.60000000e+02,
        2.16419753e+02, 1.87456790e+02, 1.70000000e+02, 1.61234568e+02,
        1.58641975e+02, 1.60000000e+02, 1.63382716e+02, 1.67160494e+02,
        1.70000000e+02, 1.70864198e+02, 1.69012346e+02, 1.64000000e+02,
        1.55679012e+02, 1.44197531e+02, 1.30000000e+02, 1.13827160e+02,
        9.67160494e+01, 8.00000000e+01, 6.53086420e+01, 5.45679012e+01,
        5.00000000e+01, 5.41234568e+01, 6.97530864e+01, 1.00000000e+02,
        1.48271605e+02, 2.18271605e+02, 3.14000000e+02],
       [5.42419753e+02, 4.30839506e+02, 3.45135802e+02, 2.81308642e+02,
        2.35654321e+02, 2.04765432e+02, 1.85530864e+02, 1.75135802e+02,
        1.71061728e+02, 1.71086420e+02, 1.73283951e+02, 1.76024691e+02,
        1.77975309e+02, 1.78098765e+02, 1.75654321e+02, 1.70197531e+02,
        1.61580247e+02, 1.49950617e+02, 1.35753086e+02, 1.19728395e+02,
        1.02913580e+02, 8.66419753e+01, 7.25432099e+01, 6.25432099e+01,
        5.88641975e+01, 6.40246914e+01, 8.08395062e+01, 1.12419753e+02,
        1.62172840e+02, 2.33802469e+02, 3.31308642e+02],
       [5.64567901e+02, 4.50617284e+02, 3.62691358e+02, 2.96790123e+02,
        2.49209877e+02, 2.16543210e+02, 1.95679012e+02, 1.83802469e+02,
        1.78395062e+02, 1.77234568e+02, 1.78395062e+02, 1.80246914e+02,
        1.81456790e+02, 1.80987654e+02, 1.78098765e+02, 1.72345679e+02,
        1.63580247e+02, 1.51950617e+02, 1.37901235e+02, 1.22172840e+02,
        1.05802469e+02, 9.01234568e+01, 7.67654321e+01, 6.76543210e+01,
        6.50123457e+01, 7.13580247e+01, 8.95061728e+01, 1.22567901e+02,
        1.73950617e+02, 2.47358025e+02, 3.46790123e+02],
       [5.80000000e+02, 4.63827160e+02, 3.73827160e+02, 3.06000000e+02,
        2.56641975e+02, 2.22345679e+02, 2.00000000e+02, 1.86790123e+02,
        1.80197531e+02, 1.78000000e+02, 1.78271605e+02, 1.79382716e+02,
        1.80000000e+02, 1.79086420e+02, 1.75901235e+02, 1.70000000e+02,
        1.61234568e+02, 1.49753086e+02, 1.36000000e+02, 1.20716049e+02,
        1.04938272e+02, 9.00000000e+01, 7.75308642e+01, 6.94567901e+01,
        6.80000000e+01, 7.56790123e+01, 9.53086420e+01, 1.30000000e+02,
        1.83160494e+02, 2.58493827e+02, 3.60000000e+02],
       [5.88567901e+02, 4.70320988e+02, 3.78395062e+02, 3.08790123e+02,
        2.57802469e+02, 2.22024691e+02, 1.98345679e+02, 1.83950617e+02,
        1.76320988e+02, 1.73234568e+02, 1.72765432e+02, 1.73283951e+02,
        1.73456790e+02, 1.72246914e+02, 1.68913580e+02, 1.63012346e+02,
        1.54395062e+02, 1.43209877e+02, 1.29901235e+02, 1.15209877e+02,
        1.00172840e+02, 8.61234568e+01, 7.46913580e+01, 6.78024691e+01,
        6.76790123e+01, 7.68395062e+01, 9.80987654e+01, 1.34567901e+02,
        1.89654321e+02, 2.67061728e+02, 3.70790123e+02],
       [5.90419753e+02, 4.70246914e+02, 3.76543210e+02, 3.05308642e+02,
        2.52839506e+02, 2.15728395e+02, 1.90864198e+02, 1.75432099e+02,
        1.66913580e+02, 1.63086420e+02, 1.62024691e+02, 1.62098765e+02,
        1.61975309e+02, 1.60617284e+02, 1.57283951e+02, 1.51530864e+02,
        1.43209877e+02, 1.32469136e+02, 1.19753086e+02, 1.05802469e+02,
        9.16543210e+01, 7.86419753e+01, 6.83950617e+01, 6.28395062e+01,
        6.41975309e+01, 7.49876543e+01, 9.80246914e+01, 1.36419753e+02,
        1.93580247e+02, 2.73209877e+02, 3.79308642e+02],
       [5.86000000e+02, 4.64049383e+02, 3.68716049e+02, 2.96000000e+02,
        2.42197531e+02, 2.03901235e+02, 1.78000000e+02, 1.61679012e+02,
        1.52419753e+02, 1.48000000e+02, 1.46493827e+02, 1.46271605e+02,
        1.46000000e+02, 1.44641975e+02, 1.41456790e+02, 1.36000000e+02,
        1.28123457e+02, 1.17975309e+02, 1.06000000e+02, 9.29382716e+01,
        7.98271605e+01, 6.80000000e+01, 5.90864198e+01, 5.50123457e+01,
        5.80000000e+01, 7.05679012e+01, 9.55308642e+01, 1.36000000e+02,
        1.95382716e+02, 2.77382716e+02, 3.86000000e+02],
       [5.76049383e+02, 4.52469136e+02, 3.55654321e+02, 2.81604938e+02,
        2.26617284e+02, 1.87283951e+02, 1.60493827e+02, 1.43432099e+02,
        1.33580247e+02, 1.28716049e+02, 1.26913580e+02, 1.26543210e+02,
        1.26271605e+02, 1.25061728e+02, 1.22172840e+02, 1.17160494e+02,
        1.09876543e+02, 1.00469136e+02, 8.93827160e+01, 7.73580247e+01,
        6.54320988e+01, 5.49382716e+01, 4.75061728e+01, 4.50617284e+01,
        4.98271605e+01, 6.43209877e+01, 9.13580247e+01, 1.34049383e+02,
        1.95802469e+02, 2.80320988e+02, 3.91604938e+02],
       [5.61604938e+02, 4.36543210e+02, 3.38395062e+02, 2.63160494e+02,
        2.07135802e+02, 1.66913580e+02, 1.39382716e+02, 1.21728395e+02,
        1.11432099e+02, 1.06271605e+02, 1.04320988e+02, 1.03950617e+02,
        1.03827160e+02, 1.02913580e+02, 1.00469136e+02, 9.60493827e+01,
        8.95061728e+01, 8.09876543e+01, 7.09382716e+01, 6.00987654e+01,
        4.95061728e+01, 4.04938272e+01, 3.46913580e+01, 3.40246914e+01,
        4.07160494e+01, 5.72839506e+01, 8.65432099e+01, 1.31604938e+02,
        1.95876543e+02, 2.83061728e+02, 3.97160494e+02],
       [5.44000000e+02, 4.17604938e+02, 3.18271605e+02, 2.42000000e+02,
        1.85086420e+02, 1.44123457e+02, 1.16000000e+02, 9.79012346e+01,
        8.73086420e+01, 8.20000000e+01, 8.00493827e+01, 7.98271605e+01,
        8.00000000e+01, 7.95308642e+01, 7.76790123e+01, 7.40000000e+01,
        6.83456790e+01, 6.08641975e+01, 5.20000000e+01, 4.24938272e+01,
        3.33827160e+01, 2.60000000e+01, 2.19753086e+01, 2.32345679e+01,
        3.20000000e+01, 5.07901235e+01, 8.24197531e+01, 1.30000000e+02,
        1.96938272e+02, 2.86938272e+02, 4.04000000e+02],
       [5.24864198e+02, 3.97283951e+02, 2.96913580e+02, 2.19753086e+02,
        1.62098765e+02, 1.20543210e+02, 9.19753086e+01, 7.35802469e+01,
        6.28395062e+01, 5.75308642e+01, 5.57283951e+01, 5.58024691e+01,
        5.64197531e+01, 5.65432099e+01, 5.54320988e+01, 5.26419753e+01,
        4.80246914e+01, 4.17283951e+01, 3.41975309e+01, 2.61728395e+01,
        1.86913580e+01, 1.30864198e+01, 1.09876543e+01, 1.43209877e+01,
        2.53086420e+01, 4.64691358e+01, 8.06172840e+01, 1.30864198e+02,
        2.00617284e+02, 2.93580247e+02, 4.13753086e+02],
       [5.06123457e+02, 3.77506173e+02, 2.76246914e+02, 1.98345679e+02,
        1.40098765e+02, 9.80987654e+01, 6.92345679e+01, 5.06913580e+01,
        3.99506173e+01, 3.47901235e+01, 3.32839506e+01, 3.38024691e+01,
        3.50123457e+01, 3.58765432e+01, 3.56543210e+01, 3.39012346e+01,
        3.04691358e+01, 2.55061728e+01, 1.94567901e+01, 1.30617284e+01,
        7.35802469e+00, 3.67901235e+00, 3.65432099e+00, 9.20987654e+00,
        2.25679012e+01, 4.62469136e+01, 8.30617284e+01, 1.36123457e+02,
        2.08839506e+02, 3.04913580e+02, 4.28345679e+02],
       [4.90000000e+02, 3.60493827e+02, 2.58493827e+02, 1.80000000e+02,
        1.21308642e+02, 7.90123457e+01, 5.00000000e+01, 3.14567901e+01,
        2.08641975e+01, 1.60000000e+01, 1.49382716e+01, 1.60493827e+01,
        1.80000000e+01, 1.97530864e+01, 2.05679012e+01, 2.00000000e+01,
        1.79012346e+01, 1.44197531e+01, 1.00000000e+01, 5.38271605e+00,
        1.60493827e+00, 0.00000000e+00, 2.19753086e+00, 1.01234568e+01,
        2.60000000e+01, 5.23456790e+01, 9.19753086e+01, 1.48000000e+02,
        2.23827160e+02, 3.23160494e+02, 4.50000000e+02],
       [4.79012346e+02, 3.48765432e+02, 2.46172840e+02, 1.67234568e+02,
        1.08246914e+02, 6.58024691e+01, 3.67901235e+01, 1.83950617e+01,
        8.09876543e+00, 3.67901235e+00, 3.20987654e+00, 5.06172840e+00,
        7.90123457e+00, 1.06913580e+01, 1.26913580e+01, 1.34567901e+01,
        1.28395062e+01, 1.09876543e+01, 8.34567901e+00, 5.65432099e+00,
        3.95061728e+00, 4.56790123e+00, 9.13580247e+00, 1.95802469e+01,
        3.81234568e+01, 6.72839506e+01, 1.09876543e+02, 1.69012346e+02,
        2.48098765e+02, 3.50839506e+02, 4.81234568e+02],
       [4.75975309e+02, 3.45135802e+02, 2.42098765e+02, 1.62864198e+02,
        1.03728395e+02, 6.12839506e+01, 3.24197531e+01, 1.43209877e+01,
        4.46913580e+00, 6.41975309e-01, 9.13580247e-01, 3.65432099e+00,
        7.53086420e+00, 1.15061728e+01, 1.48395062e+01, 1.70864198e+01,
        1.80987654e+01, 1.80246914e+01, 1.73086420e+01, 1.66913580e+01,
        1.72098765e+01, 2.01975309e+01, 2.72839506e+01, 4.03950617e+01,
        6.17530864e+01, 9.38765432e+01, 1.39580247e+02, 2.01975309e+02,
        2.84469136e+02, 3.90765432e+02, 5.24864198e+02],
       [4.84000000e+02, 3.52716049e+02, 2.49382716e+02, 1.70000000e+02,
        1.10864198e+02, 6.85679012e+01, 4.00000000e+01, 2.23456790e+01,
        1.30864198e+01, 1.00000000e+01, 1.11604938e+01, 1.49382716e+01,
        2.00000000e+01, 2.53086420e+01, 3.01234568e+01, 3.40000000e+01,
        3.67901235e+01, 3.86419753e+01, 4.00000000e+01, 4.16049383e+01,
        4.44938272e+01, 5.00000000e+01, 5.97530864e+01, 7.56790123e+01,
        1.00000000e+02, 1.35234568e+02, 1.84197531e+02, 2.50000000e+02,
        3.36049383e+02, 4.46049383e+02, 5.84000000e+02],
       [5.06493827e+02, 3.74913580e+02, 2.71432099e+02, 1.92049383e+02,
        1.33061728e+02, 9.10617284e+01, 6.29382716e+01, 4.58765432e+01,
        3.73580247e+01, 3.51604938e+01, 3.73580247e+01, 4.23209877e+01,
        4.87160494e+01, 5.55061728e+01, 6.19506173e+01, 6.76049383e+01,
        7.23209877e+01, 7.62469136e+01, 7.98271605e+01, 8.38024691e+01,
        8.92098765e+01, 9.73827160e+01, 1.09950617e+02, 1.28839506e+02,
        1.56271605e+02, 1.94765432e+02, 2.47135802e+02, 3.16493827e+02,
        4.06246914e+02, 5.20098765e+02, 6.62049383e+02],
       [5.47160494e+02, 4.15432099e+02, 3.11950617e+02, 2.32716049e+02,
        1.74024691e+02, 1.32469136e+02, 1.04938272e+02, 8.86172840e+01,
        8.09876543e+01, 7.98271605e+01, 8.32098765e+01, 8.95061728e+01,
        9.73827160e+01, 1.05802469e+02, 1.14024691e+02, 1.21604938e+02,
        1.28395062e+02, 1.34543210e+02, 1.40493827e+02, 1.46987654e+02,
        1.55061728e+02, 1.66049383e+02, 1.81580247e+02, 2.03580247e+02,
        2.34271605e+02, 2.76172840e+02, 3.32098765e+02, 4.05160494e+02,
        4.98765432e+02, 6.16617284e+02, 7.62716049e+02],
       [6.10000000e+02, 4.78271605e+02, 3.74938272e+02, 2.96000000e+02,
        2.37753086e+02, 1.96790123e+02, 1.70000000e+02, 1.54567901e+02,
        1.47975309e+02, 1.48000000e+02, 1.52716049e+02, 1.60493827e+02,
        1.70000000e+02, 1.80197531e+02, 1.90345679e+02, 2.00000000e+02,
        2.09012346e+02, 2.17530864e+02, 2.26000000e+02, 2.35160494e+02,
        2.46049383e+02, 2.60000000e+02, 2.78641975e+02, 3.03901235e+02,
        3.38000000e+02, 3.83456790e+02, 4.43086420e+02, 5.20000000e+02,
        6.17604938e+02, 7.39604938e+02, 8.90000000e+02]])

x, y = np.mgrid[-5:5:100j, -5:5:100j]
fig1,ax1=plot_surface(x,y,himmel(x,y))

np.min(np.min(himmel(x,y)))
0.0
(2.0/(1+np.sqrt(10)))**2
0.23088615702040688
(2.0/(1+np.sqrt(10)))**2
((1-np.sqrt(0.1))/(1+np.sqrt(0.1)))**2
0.26987386361223825
(2/(1+np.sqrt(0.1)))**2
2.3088615702040696
((1-np.sqrt(0.1))/(1+np.sqrt(0.1)))**2
0.26987386361223825
((np.sqrt(10)-1)/(np.sqrt(10)+1))**2
0.26987386361223836